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We describe a stochastic series expansion (SSE) quantum Monte Carlo method for a two- 
dimensional 5=1/2 XY-model (or, equivalently, hard-core bosons at half-filling) which in addition 
to the standard pair interaction J includes a four-particle term K that flips spins on a square plaque- 
tte. The model has three ordered ground state phases; for K/J < 8 it has long-range xy spin order 
(superfluid bosons), for K/J > 15 it has staggered spin order in the z direction (charge-density- 
wave), and between these phases it is in a state with columnar order in the bond and plaquette 
energy densities. We discuss an implementation of directed-loop updates for the SSE simulations 
of this model and also introduce a "multi-branch" cluster update which significantly reduces the 
autocorrelation times for large K/J. In addition to the pure J-K model, which in the z basis has 
only off-diagonal terms, we also discuss modifications of the algorithm needed when various diagonal 
interactions are included. 
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I. INTRODUCTION 

In the ongoing quest to explore possible ground states 
and quantum phase transitions in quantum condensed 
matter systems (fermions, bosons, or quantum spins), 
numerical studies are important for establishing the true 
nature of the phases and transitions of relevant model 
Hamiltonians. In particular, recent interest in "exotic" 
phenomena has focused attention on models with frus- 
trated or competing interactions, in which interplay be- 
tween adjacent ordered phases often gives rise to interest- 
ing effects j 1 ' 2 : 3 : 4 For classical models, Monte Carlo simu- 
lations in combination with finite-size scaling can be used 
very successfully in studies of a wide range of systems 
with and without frustration. However, only a limited 
class of quantum models are amenable to such studies, as 
the infamous sign problem prohibits large-scale quantum 
Monte Carlo (QMC) studies of frustrated antiferromag- 
netic spin systems and fermions in more than one dimen- 
sion. It is therefore important to search for non-sign- 
problematic quantum models, possibly with competing 
interactions, that display complex ground state phase di- 
agrams and can be efficiently studied using Monte Carlo 
simulations. Although not all possible types of ground 
states and quantum phase transitions may be realizable 
within this class of Hamiltonians, it is likely that many in- 
sights into the low-temperature physics of quantum mat- 
ter can still be gained in this way. Constructing opti- 
mized and efficient quantum Monte Carlo algorithms for 
such candidate Hamiltonians is hence an important task. 

In this paper, we present the details of a stochastic 
series expansion (SSE) algorithm that we have developed 
for large-scale QMC studies of a two-dimensional (2D) 
5 = 1/2 XY model with an added four-site ring-exchange 
term (the method can be easily generalized for three- 
dimensional systems 5 ). Defining the following bond and 



plaquette operators; 

Ba = SfSj + Srsj = 2(S?Sf + SfS'«), (1) 
Pijki = S^Sj S^S*; + Si S/j~S k S t , (2) 

the J-K Hamiltonian is given by 

H = -jJ^B ij -Kj^P m , (3) 

{ij) {ijkl) 

where {ij) denotes a pair of nearest-neighbor sites on a 
2D square lattice and {ijkl) are sites on the corners of 
a plaquette, as illustrated in Fig. da). The plaquctte- 
flip Pijki is only a subset of all the possible cyclic ex- 
changes among four spins and corresponds to retaining 
only the purely x- and y-terms; it has a non-vanishing 
matrix element only between the two spin states with 
alternating (staggered) spins on the corners of the pla- 
quette, as illustrated in Fig. ^Jb) . In the standard way, 
the J-K Hamiltonian © can also be considered as a half- 
filled hard-core boson model, where up and down spins 
correspond to filled and empty sites and J is the nearest- 
neighbor hopping. We will frequently use terminology 
referring to this boson representation. With the nega- 
tive sign in front of the plaquette-term (K > 0), the J-K 
model can be studied using QMC methods without a sign 
problem (the sign of the J-term is actually irrelevant in 
this regard). In this model, there is no frustration in the 
conventional sense, i.e., antiferromagnetic interactions on 
lattice loops with an odd number of links (which leads to 
sign problems). However, the J- and X-terms individu- 
ally favor different types of ground states, which leads to 
interesting competition effects at intermediate K/ J. 

The J-K model was recently found to exhibit three 
different ordered ground states as a function of the ra- 
tio K/J of the four-site [K) and two-site (J) terms£ It 
was argued that the transition between the magnetically 
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FIG. 1: (a) Labeling convention for the indices of an oper- 
ator Oijki acting on the corners of a plaquette. The label p 
refers to the whole plaquette, so that O p = Oijki- (b) The 
two plaquette configurations between which the if -term can 
act; open and solid circles correspond to up and down spins, 
respectively. 



updates for sampling the SSE configurations; a directed- 
loop update as well as a "multi-branch" cluster update. 
The latter significantly reduces the autocorrelation times 
for large K/J. In Sec. II we also discuss estimators for 
several important physical quantities. We discuss auto- 
correlation functions in the different ordered phases in 
Sec. IIIII In Sec. IV we discuss modifications of the algo- 
rithm when different types of potential-energy terms are 
included in addition to the J and K terms. We conclude 
with a brief discussion in Sec. V. 



ordered state for K/J < 8 and a striped (or valence- 
bond-solid, VBS) phase at higher K/J is a continuous 
quantum phase transition, contrary to general expecta- 
tions for an order-order transition. Subsequently, this 
transition was proposed to possibly be a realization of 
a "deconfined" quantum-critical point i We have used 
the SSE algorithms to further study the quantum-critical 
scaling and finite-T transitions in this model. However, 
in this paper we only briefly summarize the results and 
focus on the algorithmic issues. A full account of the 
results will be presented elsewherei 7 - 

For K = 0, the J-K model reduces to the standard 
XY-model, which undergoes a Kosterlitz-Thouless tran- 
sition at Tkt/J ~ 0.69i s i 9 i 10 In the boson language, the 
system is a superfluid below Tkt- The main features of 
the T = phase diagram for K/J > were presented 
in Ref. 6. Our most recent simulations^ 7 - show that the 
superfluid density vanishes at K/J w 7.91. At the same 
point, within the accuracy of our calculations, the ground 
state develops a stripe order, where the bond and plaque- 
tte strengths (Bij) and (Pijki) are modulated at wave- 
vector q = (it, 0) or (0, ir). This state can also be con- 
sidered a columnar VBS, since not all the bonds within 
the "ladders" of strong bonds are equal — the strongest 
ones are those on the rungs of the ladders. The VBS or- 
der vanishes at K/J w 14.5, in a first-order transition to 
an Ising-type antiferromagnetic state (a charge-density- 
wave, CDW, at q = (n, n) in the boson picture). We 
have not observed any signs of first-order behavior at the 
superfluid- VBS transition, nor any region of coexistence 
of the two phases. Numerically we can of course never 
exclude an extremely weakly first-order transition or a 
very narrow coexistence region. At the transition, we do 
observe power-law scaling with nontrivial exponents for 
the superfluid density as well as for the order parameter 
corresponding to the VBS phase. We have also recently 
studied the evolution of the VBS phase boundaries when 
the system is coupled to an external magnetic field 

The outline of the rest of this paper is the following: In 
Sec. II we describe the SSE algorithm for the J-K Hamil- 
tonian. Implementations of the SSE scheme for various 
spinels anc j -^rj f erm i on i4 models have been discussed at 
length in several recent papers, but since the four-particle 
term necessitates a more complex sampling scheme, with 
some important new features, we describe our algorithm 
in detail here. We have constructed two types of cluster 



II. STOCHASTIC SERIES EXPANSION 

The SSE method 15 i 16 i 17 i 18 is an efficient and widely 
applicable generalization of Handscomb's 1 - power-series 
method for the S = 1/2 Heisenberg model. It has 
previously been used for several models with two-body 
interactions, including the pure XY-model [K = 
in Eq. ©]- 10 As in world- line Monte Carlo, 20 loop- 
cluster algorithms^! can speed up SSE simulations very 
significantly]! 8 . Recently, a framework was devised for 
constructing and optimizing loop-type algorithms under 
very general conditions^ Here we will apply this directed 
loop scheme to SSE simulations including the four-spin 
term. A loop-type algorithm cannot be constructed for 
the pure if -model [J — in Eq. J3J], however, and the 
loops are also inefficient when J/K <C i. We therefore 
also develop a new type of multi-branch cluster update, 
that can be used in combination with the directed loops, 
and enables efficient simulations for any J/K. The multi- 
branch update bares some resemblance to, but is more 
complex than, a "quantum-cluster" update recently de- 
veloped for the transverse Ising models 

Below, we give a brief general summary of the SSE 
method. We then develop the directed loop and multi- 
branch cluster updates for the J-K Hamiltonian and dis- 
cuss the SSE estimators of several important physical 
quantities. We present some illustrative results for au- 
tocorrelation functions obtained with and without the 
multi-branch cluster update before concluding with a dis- 
cussion of the directed-loop equations for the Hamilto- 
nian with various diagonal interaction included. 



A. General SSE formalism 

To construct the SSE representation of a quantum me- 
chanical expectation value at temperature T — 1/(3; 

(A) = I T r{Ae-^}, Z = Tr{e-^}, (4) 

the Hamiltonian is first written as a sum of elementary 
interactions 

t a 
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where in a chosen basis {\a}} the operators satisfy 

H tia \a) ~ \a'}, (6) 

where \a) and \a') are both basis states. The indices t 
and a refer to the operator types (various kinetic and po- 
tential terms) and the lattice units over which the interac- 
tions are summed (bonds, plaquettes, etc.). A unit oper- 
ator Ho o = 1 is also defined. Using the Taylor expansion 
of e~P H truncated at order M, the partition function can 
then be written as^ 



Z 



EE 



(3 n {M 



M\ 



M 



where Sm denotes a sequence of operator-indices; 
Sm = [ti,ai], [t2,a,2], ■ ■ ■ , [tM, am], 



(7) 



(8) 



and n denotes the number of non-[0, 0] elements in Sm 
(i.e., the actual expansion-order of the terms). The finite 
truncation M and the use of a fill-in operator Ho,o are 
not strictly necessary^ but simplify some aspects of the 
algorithm. M can be adjusted during the equilibration 
of the simulation, so that it always exceeds the highest 
power n reached; M = ^4n max , where a suitable value for 
the factor is A w 1.25. This leads to M ~ (3N, where 
N is the system volume, and the remaining truncation 
error is completely negligible. The adjustment of M has 
been discussed in more detail in Ref. 1171 

Defining a normalized state \a(p)) obtained by acting 
on | a) = |a(0)) with the first p operators in the product 
in Eq. (J7J, 



\a(p)) ~f[H Utai \a), 



(9) 



the requirement for a non-zero contribution to Z is the 
propagation periodicity \a(M)} = \a(0)). This implies 
considerable constraints on the off-diagonal operators 
in the product, and clearly the vast majority of the 
terms are zero. In an efficient SSE method, transitions 
(a, Sm) — * (a', S' M ) satisfying detailed balance should be 
attempted only within the subset of contributing config- 
urations. Although the details of such sampling proce- 
dures to some extent depend on the model under study, 
three different classes of updates are typically used. We 
here summarize these in general terms, before turning to 
the implementation for the J-K model: 

(i) The expansion order n is changed in diagonal up- 
dates, where a fill-in unit operator is replaced by a di- 
agonal operator from the sum (J5J, and vice versa, i.e., 
Ha,o Hd, a , where the type-index d corresponds to a 
diagonal operator in the basis used. 

(ii) Off-diagonal operators cannot be added and 
removed one-by-one with the periodicity constraint 
\a(M)) = |a(0)) maintained. Local updates involving 
two simultaneously replaced operators can be used for 
this purpose. 16 However, much more efficient cluster- type 



updates, which may involve a large number of operators, 
can also be constructed. 1318 Here the general strategy is 
to find a set of operators {U, ai}, such that a new valid 
configuration can be obtained by changing only the type- 
indices ti. For the J-K Hamiltonian, we will discuss two 
such updates; directed loops and multi-branch clusters. 

(iii) A third type of update is one that affects only the 
state \a). This state, which is just one out of the whole 
cycle of propagated states \a(p)}, can change also in the 
updates (ii) involving off-diagonal operators. However, at 
high temperatures many sites will frequently have no op- 
erators acting on them. The local states at these sites will 
then not be affected by the off-diagonal updates. They 
can instead be randomly modified as they do not affect 
the weight. Such state updates can improve the statistics 
at high temperatures but are often not required for the 
sampling to be ergodic. 



B. Plaquette operators 

Turning now to the J-K model, we use the standard 
^-component basis 



\a) 



N 



), 



= ±1 



(10) 



where Sf — I /2a?, on lattices with N = L x X L y sites 
(or N plaquettes). Typically we consider square lattices, 
L x = L y , but some results for rectangular, L x ^ Ly, 
systems have also been discussed^ It is convenient to 
express all interactions in the Hamiltonian in terms 
of plaquette operators, 



(11) 





— Clijki, 


#2, a 


= (J/2)BijIki, 




= {J/2)B 3k I a , 


H^ a 


= (J/2)BkiIij, 


H5, a 


= (J/2)B u I jk , 




= KPijki, 



where 7y- and Iijki are unit operators associated with 
bonds and plaquettes, respectively, and the indexing is 
defined in Fig. Up to a constant NC, the Hamilto- 
nian is then given by a sum JSJ, where the type index 
t = 1, . . . , 6, and a is the plaquette index; a = 1, . . . , N. 
As explained above, there is also a unit operator i?o,o = 
1, which is not part of the Hamiltonian but has been 
introduced only as a fill-in element for augmenting the 
operator-index sequences of length n < M in the trun- 
cated partition function (Eq. @) to M. 



C. Diagonal update 

Because there are no diagonal operators in the original 
Hamiltonian J3J, the constant operators B\ a have been 
added in order to enable diagonal updates of the form 
[0,0] <-> [l,a] in Sm- For all elements [a p ,t p ] with t p — 
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0, 1, such substitutions can be carried out sequentially for 
p = 1, . . . , M. In the — » direction, the plaquette index 
a is chosen randomly among 1 , . . . , N. The Metropolis 
acceptance probabilities arc thcr 



,17 



P([0,0]-[l,a]) - 
P([l,a] -► [0,0]) = 



NC(3 

M-n 
M-n + 1 
NC/3 ' 



(12) 
(13) 



where P > 1 should be interpreted as probability one. If 
an attempt to remove a plaquette operator, i.e., [1, a] — » 
[0,0], is not accepted, a new plaquette index a can be 
generated at random. Note that for this model, where 
the only diagonal operators are the added constants H\ ai 
it is not necessary to keep track of the propagated states 
during the diagonal update. In general, e.g., if a diagonal 
interaction is added to the Hamiltonian J2J| , the constant 
C in Eqs. I|12|) and l|13|) should be replaced by the matrix 
element {a(p)\H ljap \a(p - 1)) = {a(p)\H ljap \a{p)) of the 
diagonal operator in the propagated state at which the 
replacement is done. 



D. Linked vertices 
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FIG. 2: Examples of vertices for the J-K model. The solid and 
open circles correspond to up and down spins, respectively, 
before (beneath the bar) and after (above the bar) an operator 
has acted, (a) is one out of 16 diagonal C-vertices, (b) is one 
out of 32 J- vertices (which flip two spins), and (c) one of the 
two K- vertices (which flip all four spins). 
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In the directed loop and multi-branch cluster updates, 
which we will discuss below, it is useful to represent the 
matrix elements in Eq. as a linked lists of "vertices" ML 
The weight of a configuration (a, Sm) can be written as 

W(a,S M ) = P n{M M ; n)l f[w(p), (14) 

P =i 

where W(p) is a vertex weight, which is simply the ma- 
trix element of the corresponding plaquette operator at 
position p in Sm', 

W(p) = (a(p)\H tp , ap \a(p-l)), (15) 

which with the operators ljll|l can take the values C, J/2, 
or K. Since the loop and cluster updates are carried 
out within sectors of fixed n (only the diagonal update 
changes n), the fill-in operators i?o,o ar e not needed in 
the linked-vertex representation. A vertex represents the 
local four-spin states on plaquette a p in the matrix el- 
ement (|15|) before and after the plaquette operator has 
acted. These eight spin states constitute the legs of the 
vertex. For the J-K model, there are three classes of 
vertices, as illustrated in Fig. [21 The constant operators 
■Hi, a correspond to C-vertices (with weight C), the bond- 
flip operators i?2,o-^5,a to J-vertices (with weight J/2) 
and the plaquette-flip operators H§ a to K-vertices (with 
weight K). An example of a linked-vertex representa- 
tion of a term with three plaquette operators is shown in 
Fig- El The links connect vertex-legs on the same site, so 
that from each leg of each vertex, one can reach the next 
or previous vertex- leg on the same site (i.e., the links are 



FIG. 3: The linked-vertex representation corresponding to 
the matrix element (a|i/6,2-ff4,i-ff4,2|a) (left), where the basis 
state \a) = j IIITTI) (right). The bidirectional links are rep- 
resented as dashed lines. The site and plaquette numbering 
for the six-site lattice is shown to the right. The numbers at 
the vertex legs indicate links across the periodic propagation 
boundary, and the corresponding spins are hence the same as 
in the state \a). The numbers here also correspond to the site 
numbering of the lattice shown to the right. 

bidirectional). In cases where there is only one operator 
acting on a given site, the corresponding "before" and 
"after" legs of the same vertex are linked to each other 
(as is the case with the legs on sites 1 and 2 in Fig. EJ . 

During the simulation, the spin state |a) and the op- 
erator list Sm are stored at all times. The linked-vertex 
representation is created after each full sweep of diago- 
nal updates. After the directed loop and multi-branch 
cluster updates have been carried out, the changes are 
mapped back into a new \a) and Sm- We will not dis- 
cuss here how these data structures are implemented and 
used in practice in a computer program. The procedures 
are completely analogous to simulations with two-body 
interactions, for which an implementation was described 
in detail in Ref. 



E. Directed loops 

In the original QMC loop algorithm, 22 spins are 
flipped along a one-dimensional closed path (the 
loop) on the space-time lattice of the discretized 
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(Trotter-decomposed) or continuous 23 path-integral 
representation^. The path is self-avoiding, and a config- 
uration can be subdivided into loops that may be flipped 
independently of each other. Allowing the path to self- 
intersect and backtrack, one can construct valid algo- 
rithms for a much larger class of models. Such gen- 
eral loop-type algorithms have been constructed both 
for continuous-time world-lines (the worm algorithm 24 ) 
and for SSE (the operator- loop algorithm 18 ). The de- 
tailed balance equations — the directed loop equations — 
that must be satisfied when constructing general self- 
intersecting and back-tracking loops were recently de- 
rived within the SSE framework, and a generalization to 
the path integral representation was also shown^ Here 
we will implement the directed-loop scheme for SSE sam- 
pling of the J and K terms. 

In an SSE operator-loop algorithm, where the loops 
constitute connected strings of operators (or vertices in 
the linked- vertex representation) ^« the building of a loop 
consists of a series of steps, in each of which a vertex 
is entered at one leg (the entrance leg) and an exit leg 
is chosen according to probabilities that depend on the 
entrance leg and the spin states at all the legs. The 
entrance to the following vertex is given by the link from 
the chosen exit leg. The spins at all vertex-legs visited 
are flipped during the loop building. 

The original starting point of the loop is chosen at 
random. Two link- discontinuities are created when the 
first pair of entrance and exit spins is flipped, i.e., the 
legs to which these are linked will be in different spin 
states (this is analogous to introducing the two sources 
in the worm algorithm 2 ^). Configurations contributing 
to Z only contain links between legs in the same spin 
states. One of the discontinuities will be propagated dur- 
ing the loop-building, whereas the other one will remain 
at the original starting point. The loop closes when the 
propagating discontinuity reaches the stationary one, so 
that they annihilate each other. A new contributing con- 
figuration has then been generated. If the path is self- 
intersecting (which is not always the casei^), the changes 
in the configuration may in effect correspond to several 
disconnected loops. 

When a vertex has been entered at a given leg, the 
probabilities for choosing one out of the possible exit 
legs have to be chosen so that detailed balance is sat- 
isfied. In general, these probabilities are not unique, and 
in most cases the most evident ones involve high probabil- 
ities for bounces, where the exit and entrance legs are the 
same and the loop building hence backtracks one stepii 3 - 
It is normally 25 desirable to minimize the probability of 
bounces. The directed loop scheme^ systematizes the 
search for valid sets of exit probabilities and enables a 
minimization of the bounce probability. To construct the 
directed loop equations for the exit probabilities, weights 
are first assigned to all possible paths through a vertex 
from a given entrance leg. The sum of all these path 
weights must equal the bare vertex weight i|15|) . i.e., the 
matrix clement before the entrance and exits spins have 
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FIG. 4: Two examples of vertex-paths that are related in the 
directed loop scheme, (a) shows a process, and its reverse, 
where a C-vertex is transformed into a J-vertex. (b) shows 
two related J<->K transformations. In the loop construction 
the spin states at the entrance and exit legs are flipped. The 
spin states shown in the vertices here are those before the flips 
have been carried out. 
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FIG. 5: A closed set of C and J vertex-paths, with their 
corresponding weights ay that have to satisfy the directed 
loop equations. 



been flipped. The actual normalized exit probability is 
the path weight divided by the bare vertex weight. The 
key element of the scheme is that weights for vertex-paths 
that constitute each other's reverses have to be equal in 
order for detailed balance to be fulfilled (a generalized 
scheme where this is not necessarily the case has also been 
discussed recently 5 .). Examples of such related vertex- 
paths in the J-K model are shown in Fig. The directed 
loop equations written down on the basis of these simple 
rules often form several different closed sets that can be 
solved for the path weights independently of each other. 
Because of symmetries, many of the equation sets can 
also be also identical. In general, the directed loop equa- 
tions have an infinite number of solutions, which can be 
significantly restricted by minimizing the bounce proba- 
bilities. In some cases there is a unique minimum-bounce 
solution (sometimes with zero bounce probability), but 
often there is still a high degree of freedom \q{xM^^ 

For the J-K model, a one dimensional path segment 
can in one step transform a C-vertex into a J-vertex, and 
vice versa, an example of which is shown in Fig. Ufa). 
A J-vertex can be transformed into a K-vertex, and vice 
versa, as shown in Fig. Efb). C- and K- vertices cannot 
be directly transformed into each other, however. As a 
consequence, the closed sets of vertex-paths that contain 
C<-^J transformations are independent from those con- 
taining J<-»K transformations. 

The closed sets containing C^J transformations are 
similar to those for the XY-modelji 2 although the sets 
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are larger because a C-vertex can be transformed into 
two different J-vertices. As in the XY-model, no bounces 
are required for detailed balance in this case, until we 
discuss the inclusion of additional diagonal interactions 
in Sec. IV. One closed set with C^->J transformations is 
shown in Fig. 03 To construct such a set, one first selects 
a "reference" vertex (any vertex) and an entrance leg, 
and then finds all paths that lead to new valid vertices, 
sampling all allowed exit legs. This corresponds to the 
first row of Fig. [JJ where the bounce process has not been 
included since, as will be shown below, its weight can be 
set to zero in this case. Each of the resulting vertices (i.e., 
when the entrance and exit spins have been flipped) are 
then considered in turn, using as the entrance legs the 
exit legs from the previous step. This leads to rows two 
to four in Fig. The procedure is repeated for each 
new combination of vertex and exit leg that is created. 
This systematically generates all pairs of vertex-paths 
that constitute each other's reverses, i.e., those that must 
have equal weights for detailed balance to be satisfied. In 
the case considered here, no new vertex-paths are created 
after row four, as the reverse of each path has then al- 
ready been generated. The set is hence closed. Other 
closed sets are constructed by picking a starting vertex 
and entrance leg combination that has not yet appeared 
within the sets already completed. This is repeated un- 
til all vertices and entrance legs combinations have been 
exhausted. 

The directed loop equations corresponding to the 
closed set shown in Fig. are 



an + a 12 + a 13 
«21 + a 22 + «23 
031 + &32 + 033 
041 + 042 + 043 



Wi = c, 
w 2 = c, 

W 3 = J/2, 
Wa = J/2, 



(16) 



where the weights aij are identified with the paths in 
the figure and Wi are the bare vertex weights before the 
entrance and exit spins have been flipped. Detailed bal- 
ance requires that the weights corresponding to opposite 
vertex-paths are equal, i.e., 



It can be expected that it is efficient to maximize the 
weights of the paths that transform a C-vertex into a 
J-vertex, which is equivalent to minimizing the weights 
of the continue- straight paths that transform a C-vertex 
into another C-vertex. We have no proof of our assertion 
that this is a good strategy, but as it is a quite challenging 
task to investigate all possible valid solutions, we will 
use it and leave other possibilities for future studies (this 
issue has in fact recently been addressed in the context 
of other models^). In Fig. [5J there are only two C^C 
paths; the pair with weights an, 021- The minimum 
value of these is an = 0-21 = C — J/2, which also implies 
C > J/2. There are now enough conditions to render a 
unique solution to this set of directed loop equations; 



an 

021 
031 
041 



C - J/2, 

J/4, 

J/4, 



Ol2 


= J/4, 


O13 


= J/4, 


022 


= J/4, 


023 


= J/4, 


032 


= J/4, 


033 


= 0, 


042 


= J/4, 


043 


= 0. 



(18) 



The actual exit probabilities P?- = a^/Wi are 



pa 
r U 
pa 

21 
pa 

21 
pa 
-*41 



1-J/2C, Pf 2 
1-J/2C, P? 2 
1/2, 
1/2, 



pa 

pa 
^42 



J/ AC, Pf 3 

J/*C, P 2 a 3 

1/2, 
1/2, 



pa 
^23 
pa 

M3 



J/4C, 
J/4C, 
0, 

0, 



(19) 



where the superscript a is used as a reminder that these 
probabilities correspond to the paths shown in Fig. |S] 
Note that the probabilities here depend only on the type 
of vertex transformation, C^C (P = 1 — J/2C), C— >J 
(P = J/2C), J^J (P = 0), or J^C (P = 1/2), which 
can aid the implementation of the probability tables in 
the code. All other sets with C^J transformations are ei- 
ther related by trivial symmetries to that shown in Fig. [S] 
or are very similar to it. The exit probabilities are given 
simply by the type of the corresponding vertex transfor- 
mation exactly as above. 

The directed loop equations for the closed sets of 
paths that involve J<->K transformations sometimes re- 
quire non-zero bounce probabilities. A closed set of paths 
is shown in Fig. [3] The corresponding equations for the 
path weights bij 



(20) 



021 


= an, 
















031 


= 012, 




Oil - 


h 612 - 


1- O13 " 


1- 614 - 


F-615 


- J/2, 


032 


= 2 2, 




621 - 


h 622 - 


1- 23 " 


F" 024 ' 


F-&25 


- J/2, 


«41 


= 013, 


(17) 


631 - 


h 632 ~ 


F- 033 - 


F- 034 - 


F-&35 


= K, 


O42 


= 23 , 




641 - 


h 642 ~ 


H643- 


F- 044 - 


F- 645 


= J/2, 


O43 


= a 33 . 




051- 


h& 52 - 


K653- 


F- 054 - 


F- 055 


= J/2. 



The weights also have to be positive definite, since they 
are related to probabilities by dividing with the positive 
matrix elements Wi- Even with these constraints, the 
solution is not unique. One can reasonably assume that 
the most efficient solution also has equal weights for paths 
that are related by symmetries, e.g., a 12 = 013. Using all 
such symmetries, the solution is still not unique, however. 



Again, it is in general advantageous to minimize the 
bounce probabilities, i.e., the bounce weights 6^5 above. 
For K < 2J all the bounce weights can in fact be zero. 
The weight of the continue-straight paths (e.g., 6n), 
which here transform a J-vertex into a J-vertex with the 
same spin flips (i.e., the same plaquette operator), can 
be set to zero. A symmetric K<2J solution is then: 
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FIG. 6: A closed set of vertex-paths with J^K transformations, labeled by their path weights bi 



fell =0, b 12 = if/4, 613 = J/4 -if/8, 6 14 = J/4-if/8, b 15 = 0, 

&21 -0, 6 22 = if/4, & 23 = J/4 -if/8, & 24 = J/4 -if/8, & 25 =0, 

6 31 =if/4, 6 32 =if/4, & 33 = -fiT/4, 6 34 =if/4, & 35 -0, (21) 

641 = J/4 -if/8 fo 42 =0, 643 = J/4 -if/8, fe 44 =if/4, fe 45 =0, 

651 = J/4 -if/8, 652 =0 653 =if/4, fe 54 =J/4-if/8, 655 =0. 

For if > 2J, the bounce weight 635 has to be non-zero for a positive-definite solution. Minimizing this weight one 
obtains the following solution: 



(22) 



hi 


= 0, 


bi2 


= J/2, 


bi3 


= 0, 


&14 


= 0, 


bi 5 


= 0, 


621 


= 0, 


1)22 


= J/2, 


b 2 3 


= 0, 


^24 


= 0, 


hs 


= 0, 


631 


= J/2, 


b 3 2 


= J/2, 


&33 


= J/2, 


kii 


= J/2, 


hb 


= if - 2 J 


&41 


= 


bi2 


= 0, 


bi3 


= 0, 


644 


= J/2, 


bis 


= 0, 


&51 


= 0, 


b 5 2 


= 


&53 


= J/2, 


b 5 4 


= 0, 


&55 


= 0. 



The exit probabilities are hence, for if < 2 J: 



pb 


= 0, 




pb 

r 12 


= if/2 J, 


pb 

r 13 


= 1/2 - if /4 J, 


pb 

•VL4 


= 1/2 - if/4J, 


pb 


= 0, 


pb 

r 21 


= 0, 




pb 
r 22 


= if/2 J, 


pb 

r 23 


= 1/2 - if /4 J, 


pb 
r 2i 


= 1/2 - if/4J, 


pb 

r 25 


= 0, 


pb 
r 31 


= 1/4, 




pb 
r 32 


= 1/4, 


pb 
r 33 


= 1/4, 


pb 
r 3i 


= 1/4, 


pb 

^35 


= 0, 


pb 


= 1/2- 


K/4J 


pb 

r i2 


= 0, 


pb 

r A3 


= 1/2 - if /4 J, 


pb 


= K/2J, 


pb 


= 0, 


pb 


= 1/2- 


K/4J, 


pb 

r 52 


= 


pb 
^53 


= if/2 J, 


pb 


= 1/2 - if/4J, 


pb 


= 0, 



and for if > 2 J: 



pb 


= 0, 


pb 

•VL2 


= 1, 


P b 

r 13 


= 0, 


P b 

■V14 


= 0, 


pb 


= 0, 


pb 

r 2\ 


= 0, 


pb 

r 22 


= 1, 


pb 
r 23 


= 0, 


pb 
r 2i 


= 0, 


pb 
r 2Z 


= 0, 


pb 
^31 


= J/2if, 


pb 
r 32 


= J/2if, 


pb 
■M33 


= J/2if, 


pb 
^34 


= J/2if , 


pb 
*35 


= 1 


pb 
Ml 


= 


pb 


= 0, 


pb 


= 0, 


pb 
M4 


= 1, 


pb 
M5 


= 0, 


06 
^51 


= 0, 


pb 
r b2 


= 


pb 
^53 


= 1, 


pb 

^54 


= 0, 


pb 


= 0. 



(23) 



2J/if, (24) 



Note that the solution is continuous across if = 2 J. 



Also in this case the probabilities are seen to depend 
only on the type of vertex class transformation, J— >J, 
J— »J', J^K, K^J, or K^K (bounce). Here one has 
to distinguish between a continue-straight J— >J transfor- 
mation where the spin-flip remains on the same bond 
(e.g., &n), and a J— >J' transformations where the spin- 



flip moves to a neighboring bond on the plaquette (e.g., 
613). 

There is one more type of closed set of vertex-paths, an 
example of which is shown in Fig.0 In this case, neither 
a valid K-vertex nor a J-vertex with the flip moved to a 
different nearest-neighbor pair can be reached from the 
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FIG. 7: A closed set of two vertex-paths, where only the 
continue-straight process is allowed and a J-vertex is trans- 
formed into another J-vertex with probability 1. 

TABLE I: All exit probabilities for the J-K model. The initial 
vertex class is indicated in front of the square bracket, and 
the new class after the arrow. All the possible vertex classes 
that can be generated from a given vertex and entrance leg 
are listed within the square bracket. In cases where more than 
one vertex of a given class can be generated, the corresponding 
symbol appears multiple times. J, J', J" denote subclasses 
of J-vertices in which different spin pairs are flipped. The sets 
a,b,c correspond to Figs. |U H3 all other sets are related to 
these by symmetries. The only bounce process is K— >K; all 
C^C and J-^J cases correspond to continue-straight paths. 

Vertex transformation set P(K < 2 J) P{K > 2 J) 



C-[C,J,J']-> C a 1 - J/2C 1 - J/2C 

C-[C,J,J']-> J,J' a J/4C J/4C 

J-[C,C,J']-» C,C a 1/2 1/2 

J-[C,C,J'H J' a 

J-[J,J',J",K]-» J b 

J-[J,J',J",K]-> J', J" b 1/2- K/AJ 

J-[J,J',J",K]-> K b K/2J 1 

K-[J,J,J',J',K]-> J,J,J',J' b 1/4 J/2K 

K-[J,J,J',J',K]-> K b 1 - 2.J/K 

J-[J]-> J c 1 1 



J-vertex and the chosen entrance leg. As the two vertices 
shown have the same bare weights, no bounce processes 
have to be included and the exit is unique: 

n = i, 

P 2 C =1. (25) 

All closed sets of vertex-paths can be related to those 
shown in Figs. El and [7| and in all cases the proba- 
bilities depend only on how the paths transform the ver- 
tices between the classes C, J, and K. This simplifying 
property will be discussed further in Sec. IV, where we 
consider inclusions of additional diagonal interactions in 
the Hamiltonian. In that case the solutions to the equa- 
tions become more complicated, however the directed- 
loop framework is still required in order to develop effi- 
cient codes. For the case of zero diagonal interactions, 
Eq. J2J), the exit probabilities for the J-K model are sum- 
marized in Table 

To carry out a directed loop update, a vertex-leg is 
first chosen at random. This entrance leg together with 
all the spin states on the vertex determine which one of 
the exit probabilities in Table [I] should be applied when 
generating the exit leg. These probabilities can be stored 



in a pre-generated table. When the exit has been se- 
lected, the link from it is used to enter another vertex, 
from which an exit is again chosen, etc., until the loop 
closes. The number of loops to be generated during each 
Monte Carlo step is adjusted such that the total number 
of vertices visited is, on average, of the same order as 
(e.g., equal to or twice) the number of vertex- legs (8n), 
e.g., 4(n) or AM. 

In some cases, a loop can become very long before it 
closes. In order to avoid problems with loops that do 
not close within a reasonable time, one can impose a 
maximum loop length. If this limit is exceeded, the loop 
building is terminated and the changes in the vertices 
are disregarded. This does not introduce any bias in 
quantities measured in the (a, Sm) representation. In 
practice, the termination can easily be accomplished by 
simply exiting the loop-update routine without mapping 
the linked- vertex representation back into a state \a) and 
an operator list Sm', in order to discard only the loop 
currently under construction, its history would have to 
be stored. Hence, not only the terminated loop itself is 
discarded, but also all other loops constructed since the 
previous diagonal update. This is not a problem as long 
termination does not occur frequently. We typically set 
the maximum loop length to w 100(n), and the fraction 
of terminated loops is then very small. 



F. Multi-branch clusters 

Since a K-vertex cannot be generated directly out of 
a C-vertex, but requires the presence of J-vertices, the 
directed loop update cannot be used when J = 0. As 
will be demonstrated in Sec. IIIII it is also inefficient for 
large K / J. This can be understood from Table H where 
the bounce probability off a K-vertex is seen to approach 
1 as K — > oo. In principle the directed-loop update, 
in combination with the diagonal update, is ergodic for 
any finite K/J, but for K/J > 12 it becomes difficult 
to obtain good results this way. In order to improve the 
performance for large K/J, a type of multi-branch clus- 
ter update is developed here. It is similar to a quantum- 
cluster update recently developed for the transverse Ising 
modelfl^ where it can be considered a direct generaliza- 
tion of the classical Swendsen-Wang algorithm. 29 The 
multi-branch cluster update for the J-K model is more 
complex, due to the larger number of different interac- 
tion vertices and the multitude of possible transforma- 
tions among them. 

In order to transform a C-vertex directly into a K- 
vertex, spins at four legs have to be flipped. If this is 
done, spins also have to be flipped at all the legs to which 
these four legs are linked. This will in turn force addi- 
tional spin flips in the vertices to which they are linked, 
etc. Clearly, such a process can branch out very quickly 
to a large number of vertices. Even if a scheme can be 
found where detailed balance is maintained, there is in 
general nothing that guarantees that the process ever ter- 
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FIG. 9: Multi-branch cluster update in which two C-vertices 
are transformed into two K- vertices. The initial entrance leg 
is at the inward pointing arrow in the linked- vertex represen- 
tation to the left. The resulting vertices with their arrows 
indicating legs visited are shown to the right. 
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FIG. 8: Vertex transformations in the multi-branch cluster 
update. The entrance leg is denoted by an arrow pointing 
into the vertices on the left. In the updated vertices to the 
right, the spins at the outgoing arrows have been flipped. The 
branching for all entrance legs and vertices not shown here are 
obtained by applying trivial symmetries to one of the cases 
shown in (a)-(i). 



minates (i.e., the cluster may not complete). For the 
J-K model, this proliferation problem can be solved by 
defining a unique set of exit legs, given an entrance leg 
and all the eight spin states of the vertex. If the con- 
stant C is chosen equal to K, which with the directed 
loop probabilities in Table |U can be done if K > J/2, 
C<->K transformations will lead to no weight changes. If 
J-vertices are always transformed into other J-vertices, 
there are also no weight changes in these processes. A 
constructed cluster can therefore be flipped with proba- 
bility one. One can also subdivide the whole linked ver- 
tex list into clusters that can be flipped independently of 
each other with probability 1/2. This Swendsen- Wang- 
type approach will be used here. 

Fig. [5] shows the branching rules for all different types 
of vertices. The cases (d) and (e) correspond to C^K 
transformations. In all other cases the vertex class does 
not change, but note that J-vertices are transformed into 



J-vertices with a different pair of flipped spins (i.e., the 
corresponding plaquette operator changes). The outgo- 
ing arrows point to entrances to other vertices, to which 
the same branching rules are applied. However, if an exit 
leg is linked to a leg which has already been visited, this 
leg should not be visited again. In terms of the graphical 
representation used in Fig. [8] a vertex-leg should not be 
entered if it already has an outgoing arrow. If a vertex is 
entered for a second time, and hence has arrows at four 
legs (those with eight exit legs in Fig. [H] can clearly only 
be visited once) , the second set of exit legs are exactly the 
four that were not previously assigned arrows. In other 
words, all vertices in (c)-(i) can be assigned outgoing ar- 
rows in two different ways, and the set chosen is the one 
to which the entrance leg belongs. Furthermore, the two 
sets of mutually exclusive exit legs are exactly the same 
in the vertices obtained when the legs in one of these sets 
are flipped. This solves the proliferation problem, since it 
is guaranteed that a vertex-leg can be visited only once. 
It also allows for independent flips of all clusters. 

To start a cluster, a vertex-leg which does not belong 
to a cluster already constructed is first chosen at random, 
and the branching is assigned according to the rules de- 
fined in Fig. |S1 Flags are set on all the exit legs, to in- 
dicate that they have been visited (corresponding to the 
outgoing arrows Fig. [HJ . Note that the entrance also be- 
comes an exit leg with an outgoing arrow. If the cluster 
is to be flipped (which it should with probability 1/2), 
the spins at all the exit legs are flipped. All exit legs are 
put on a stack. They are subsequently picked one-by-one 
from the stack, and the legs to which they are linked are 
used as entrance legs to other vertices if they have not 
yet been visited, i.e., these legs are flipped and put on 
the stack only if they have not been visited before. In 
the graphical representation, a cluster-branch ends when 
an arrow is encountered. The whole cluster is completed 
when all arrows point to other arrows; the stack with 
unprocessed entrance legs is then empty. A completed 
cluster with only two vertices is illustrated in Fig. [5] 

Although the autocorrelation measurements discussed 
in Sec. IIHI provide a quantification of the "effectiveness" 
of the multi-branch cluster updates, we pause here to sim- 
ply illustrate the cluster characteristics as implemented 
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Cluster size (no. of legs) 

FIG. 10: Probability histogram of cluster sizes (measured as 
the number of vertex legs in a cluster) produced by multi- 
branch cluster updates of a 16x16 lattice at K/J — 80 and 
f3 = 32. The average length of the operator list for this sim- 
ulation was (n) « 4.35 x 10 5 (i.e., approximately 3.5 x 10 6 
vertex legs). Data in the upper figure is for the smaller bins, 
while data in the lower figure was re-binned to 100 legs per 
bin. All intermediate occupations were measured as zero. 



for our J-K model. A histogram of cluster sizes gener- 
ated at large K/J is shown in Fig.EH Clearly, the vast 
majority of clusters built in this case are small eight- 
leg clusters (an example of which is illustrated in Fig.|3J, 
with significant occupation in the smaller bins up to clus- 
ters of size 128. A second peak occurs in the histogram 
at much larger bin sizes, approaching the total num- 
ber of operators. Clearly, the efficiency of the algorithm 
would be much greater if clusters occurred with sizes that 
were more evenly distributed between 8 and 8n, but in 
any case the multi-branch update does significantly im- 
prove the performance for large K/J (as will be shown 
in Sec. III). This can probably be explained by the fact 
that the directed-loop updates become very inefficient as 
K/J — > oo, and hence the multi-branch clusters help sig- 
nificantly even though the cluster-size distribution is not 
optimal. 



G. Physical observables 

In this section we summarize the physical observables 
relevant to studies of the J-K model^SiLii and present 
the estimators used to evaluate them in the SSE method. 
The general forms of the estimators have been derived in 



previous papers ^*iSiiL^S here we only apply those de- 
rived forms to the particular quantities of interest for the 
J-K model. 

We typically carry out measurements on the configu- 
rations generated after every Monte Carlo step (MCS), 
with an MCS defined as a sweep of diagonal updates, fol- 
lowed by construction of the linked vertex list, in which 
a fixed number of loop updates are carried out. In the 
same linked list, all multi-branch clusters are constructed 
and flipped with probability 1/2. After this, the updated 
vertex list is mapped back into a new state \a) and an 
operator list Sm- This is the representation used for the 
measurements. The fill-in elements i?o,o m Sm are irrel- 
evant at this stage, and we therefore now consider the 
reduced list S n without these operators. There are hence 
n + 1 propagated states \a(p)) = \af(p),...,a^{p)), 
which are obtained one-by-one when operating with the 
first p operators, p = 0, . . . , n, on the initially stored state 
|cv(0)) = \a(n)). Although measurements can involve all 
the states, at any given time only a single |cv(p)) has to 
be stored. 

The ^-component of the spin-spin correlation function 
can be easily obtained, as it is diagonal in the representa- 
tion used. Equal-time correlations can be averaged over 
the propagated states, i.e., 

( 5 ^)=i<^£^ K(p) /' (26) 

\ P=0 / 

where in the special case n — 0, which occurs in practice 
only for small N at very high temperatures, the aver- 
aged sum should be replaced by erf (0)of (0). Since states 
p and p + 1 differ only by two or four flipped spins, the 
sum in (|26|) can be replaced by a sum where only, e.g., ev- 
ery A^ th state is included. We often consider the Fourier 
transform of the correlation function, i.e., the static spin 
structure factor 

Ss(q x ,q y ) = ^Xy ( '*~ r °' q < 5 * S *>' (27) 
k,l 

where r, = (xi,yi) is the lattice coordinate (with lattice 
spacing 1) and q = (q x x, q y ) is the wave-vector. We also 
study the corresponding static susceptibility, 



Xs{q x ,q y ) = ^E ei(rfc_r,) ' q / (Sk(r)Sf(0)). (28) 

k,l o 

It has been showniSi that the SSE estimator for the Kubo 
integral is 

/^(,)5f(0))^(^L_x (29) 
o 

(n-1 \ /n-1 \ n-1 ] \ 

E °*(p) E + E *fc(pK(p) ) • 
p=0 ) \p=0 ) p=0 J / 
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Here the first term typically dominates; it is obtained by 
first summing the spins at k and I over the propagated 
states, and then multiplying the sums. The full sums 
must clearly be calculated here, but one can still take 
advantage of the fact that only two or four out of the N 
spins a%, (p) change at every propagation p — ■> p + 1 . One 
can thus evaluate the sums for all sites k in ~ n ~ N/3 
steps. The second term in <|3(Jll vanishes as N — > oo, 
but typically it gives a non-negligible relative contribu- 
tion for small N calculations and should always be kept. 
This sum is the same as in the equal-time correlation 
(|26[) and can again be replaced by a partial summation 
without introducing a bias. In the case n = 0, the whole 
expression within (} in Eq. l|3"U)) should be replaced by 

We are also interested in the spin stiffness, or the su- 
perfluid density in the boson representation, which at 
T = is defined by 



d 2 E(cj)) 



(30) 



where E(<j>) = (H(4>))/L 2 is the ground state energy per 
site and <j> is a twist which is imposed on all bonds 
in either the x or y lattice direction, so that the corre- 
sponding bond operators become 



+ S in(<P)(SfS) 



S?SV) 
- S?S*) 



(31) 



This leads to a shift in the ground state energy E to sec- 
ond order in <fr. With the plaquette operator Pijki(4>), 
the leading-order energy shift is oc 4> A , and hence it will 
not appear in the estimator for the stiffness. The deriva- 
tive at 4> = in Eq. (|30l) can therefore be directly esti- 
mated using the winding number fluctuations in the SSE 
simulations^ in a way very similar to the way it is done 
in path integral methods^ Defining the winding num- 



bers w x and w y as 



(N+ ~ N+)/L, 



(32) 



where denote the number of operators in S n which 
transport a boson (or spin-f) ±1 lattice steps in the in- 
direction. In the J-K model, only the bond operators Bij 
can transfer a net number of particles; in terms of the 
corresponding plaquette operators l|llfl . the pairs i?2,a, 
H^ a and H 3 a , H 5 a transfer particles along the x- and 
y-axis, respectively. By operating successively with all 
operators in S n on the state \a) one can determine all 
the numbers needed to obtain the winding numbers. 
The stiffness is then given by 



2/3 



(w 2 x +wl). 



(33) 



At finite T, the ground state energy E(<f>) in Eq. I|3(J|) 
should be replaced by the free energy F(4>). It turns out 
that this leads to exactly the same estimator, Eq. I|33|) . A 



detailed derivation of this well known result for lattice 
models has been presented in Ref. 13^ . 

In order to detect the modulations of the bond and 
plaquette expectation values (Bij) and (Pijki) in the 
striped phase, one can use open boundary conditions in 
order to break the translational symmetry. In order to 
break the 90° rotational symmetry, rectangular lattices 
can be used. On these lattices one can observe a unique 
bond/plaquette pattern^ However, for careful finite-size 
scaling studies it is preferable to consider periodic L x L 
lattices, on which all bond and plaquette expectations 
average to uniform values. We hence instead consider 
the corresponding correlation functions, and also calcu- 
late the associated susceptibilities. The static plaquette 
structure factor is defined as 



^(fe,^)-^E el(ra " rb) ' q ( p ^)' 

a, b 



(34) 



where P a is the plaquette operator @ with the plaquette 
subscript a defined in Fig. ^ The corresponding suscep- 
tibility is completely analogous to Eq. J33J), 



X P (Qx,q y ) 



N ^ 



dT(P a (T)P b {0)). (35) 



Bond structure factors and susceptibilities are defined in 
the same way; we here consider those corresponding to 
correlations between bonds in the same lattice direction. 
Hence, defining Xk and yk as the nearest-neighbor sites 
of site k in the x- and y-directions, the bond structure 
factors Sb, x and Sf, t y are 

S b ,*{q x ,q v ) = ^Y. STk ' n) " i ^ h B l , m ), (36) 

k,i 

and clearly S b , x (q x ,q y ) = S^ y (q y ,q x ). The corresponding 
susceptibilities are again defined as in Eq. (|35|) . 

For expectation values involving products of operators 
that also appear as terms in the Hamiltonian, such as 
the above plaquette and bond structure factors and sus- 
ceptibilities, the SSE estimators are remarkably simple 
expressions involving only numbers of operators or oper- 
ator combinations in the list S n i& The simplest case is 
the expectation value of a single operator, 



(Ht,a) 



(37) 



where n([a, b]) is the number of elements [a, b] in the list 
S n . This gives the internal energy 



E = 



(38) 



which is identical to the expression obtained by 
Handscomb^S An equal-time correlation function of two 
operators appearing in the Hamiltonian is given byA& 

(H s<a H t , b ) = jp((n- l)N([s, a][t, b})), (39) 
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FIG. 11: Autocorrelation function for the spin stiffness (upper 
panel) and the plaquette-stripe order parameter (lower panel) 
in simulations of L = 16 systems at Kj J = 4 and 7, both at 
inverse temperature K/T = 16. Results of simulations both 
with and without the multi-branch cluster update are shown. 



where N([s, a][t,b]) denotes the number of occurrences 
of the operators [s, a] and [t, b] next to each other, in the 
given order, in S n (with the periodicity of S n taken into 
account). The corresponding Kubo integral is— 
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FIG. 12: Autocorrelation function for the plaquette-stripe or- 
der parameter at K/J = 12 and inverse temperature K/T = 
32. Results with (solid curves) and without (dashed curves) 
multi-branch clusters are compared for three different system 



dT(H s>a (r)H t J0)) = (40) 



-(N([s,a])N([t,b}) - S st 5 ab N([s,a})), 
P 

where iV([s,a]) is the number of operators [s,a\. Using 
Eqs. iJSSl and JUJ, the estimators for pi jl -lp ^l can be 
easily obtained. 



III. AUTOCORRELATIONS 

We here show some results illustrating the performance 
of the algorithm, focusing in particular on the efficiency 
boost achieved with the multi-branch update. It would 
clearly be interesting to extract the dynamic exponent 
of the simulations at the various phase transitions, but 
we will not attempt this here. Instead, we will focus 
on the simulation dynamics inside the ordered phases. 
Particularly in the striped and staggered phases, which 
break spatial symmetries, we expect slow modes corre- 
sponding to transitions between the different degenerate 



states. One might also expect potential problems related 
to long-lived defects forming in these states. 

For a quantity Q, the normalized autocorrelation func- 
tion is defined in the standard way as 



A[Q](t) 



(Q(i + t)Q(i))-(Q(i))< 

(QW 2 ) 



(41) 



where the averages are over the Monte Carlo time (steps) 
i. We will compare autocorrelation functions in the three 
different ordered phases, obtained in simulations with 
and without multi-branch cluster updates. A Monte 
Carlo step is defined as a full sweep of diagonal updates, 
followed by a number of directed-loop updates, and, if 
multi-branch updates are carried out, decomposition of 
the configuration into clusters, each of which is flipped 
with probability 1/2. In these simulations the number 
of directed-loop updates per step was chosen so that, on 
average, the total number of vertices visited is 4M, with 
the truncation M of the index sequence chosen equal to 
1.25 times the maximum expansion order n reached dur- 
ing equilibration (the dependence of M on the length of 
the equilibration is in practice very small and introduces 
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FIG. 13: Autocorrelation function for the staggered order 
parameter at K/J = 32 and inverse temperature K/T = 32. 
Results with (solid curves) and without (dashed curves) multi- 
branch clusters are compared for three different system sizes. 



only a negligible ambiguity in the definition of the Monte 
Carlo time). 

Fig. ^] shows autocorrelation results for the super- 
fluid density p s and the squared stripe-order-parameter 
Mp inside the superfluid phase for a 16 x 16 lattice at 
K/T = 16. At K/J = 4, the p s autocorrelations drop 
very rapidly (the integrated autocorrelation time is less 
than 1), and there are no discernible effects of includ- 
ing multi-branch updates. The autocorrelation time for 
Mp is also very short, but here there are clear improve- 
ments with the multi-branch updates. However, consid- 
ering that the CPU time is almost doubled when includ- 
ing multi-branch updates, including them at K/J = 4 is 
not advantageous. At Kj J — 7, which is approaching 
the transition point to the striped phase at K/J w 7.9, 
the autocorrelations decay much slower, and although 
there are visible favorable effects of the multi-branch up- 
dates in both quantities, the gain is hardly worth the 
additional CPU time cost. In Fig. results are shown 
for the stripe order parameter at K/J = 12, well in- 
side the striped phase, for three different system sizes 
at inverse temperature K/T — 32. The multi-branch 
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FIG. 14: The staggered magnetization vs the propagation 
index p (divided by the total number of operators n) in con- 
figurations generated for system sizes L = 8, 16, and 32, at 
K/J = 32, K/T = 32. 



updates have clear favorable effects on the autocorrela- 
tions, but although the initial drop is considerably faster, 
the asymptotic autocorrelation time, i.e., the long-time 
linear decay seen on the linear-log scale used in the fig- 
ures, changes very little. In this case the reduction of 
the integrated autocorrelation time may (depending on 
the exact value of K/J, and the system size) motivate 
the additional computational effort of the multi-branch 
update. 

The multi-branch cluster update improves the simu- 
lation efficiency considerably inside the CDW phase, as 
illustrated in Fig. ^| for three different system sizes with 
K/J = 32 at a low temperature. Here the improvement 
in the simulation efficiency for the squared staggered or- 
der parameter M| is clearly significant enough to moti- 
vate the cost of the multi-branch clusters, especially for 
large system sizes. An interesting feature to note here 
is that when the multi-branch clusters are included, the 
asymptotic autocorrelation time actually decreases for 
L = 32 relative to L = 16, and L = 16 and L = 8 show 
almost identical autocorrelation functions. This surpris- 
ing trend for increasing L can probably be traced to the 
fluctuations in the CDW order parameter for a given SSE 
configuration. Fig. 1141 shows the dependence of the stag- 
gered order parameter on the propagation number p [re- 





time [10 Monte Carlo steps] 

FIG. 16: Bin averages for the superfluid stiffness, the stripe 
structure factor, and the stripe susceptibility, for an L — 96 
lattice at Kj J = 7.91 and inverse temperature K/T — 64. 
Each point represents an average over 10 4 Monte Carlo steps. 



FIG. 15: Distribution of the staggered magnetization at 
K/J = 32, K/T = 32. 



fcrring to the propagated states, Eq. JJJJ] divided by the 
total number of operators n for an equilibrated configu- 
ration. The fraction p/n corresponds roughly to the nor- 
malized imaginary-time t//3 in the standard Euclidean 
path integral formalism^ For a small system, exempli- 
fied by L = 8 in the figure, the order parameter fluctu- 
ates between positive and negative values, whereas for 
a large system, exemplified here by L = 32, fluctuations 
sufficiently large to "tunnel" the system between positive 
and negative order parameters are very rare. Clearly, as 
T — > 0, there would be such tunneling events also in a 
large system, but if T is not low compared to the gap 
between the symmetric and antisymmetric linear com- 
binations of the two different real-space ordered states 
(which decreases exponentially fast with increasing L), 
such events are not present in typical configurations. The 
shorter autocorrelation time for L = 32 than for L = 16 
(when multi-branch updates are included) in Fig. lI3l mav 
hence be related to the larger fluctuations in Ms for the 
smaller system size, which can lead to various tunnel- 
ing events that are not so easily added or removed from 
the configurations. For the larger system size, there are 
in practice no tunneling events at the temperature used 
here, and the difficulties in adding/removing them in this 
case would only show up at very long times as an unde- 



tectable tail in the autocorrelation function. The order 
parameter fluctuations are further illustrated in the form 
of histograms in Fig. 1 1 51 (only showing the positive part 
of the distributions). Here it can be seen that some tun- 
neling events are still expected at L = 16, since the prob- 
ability of a zero order parameter is not negligible, but at 
L = 32 the probability at zero is exponentially small and 
hardly any tunneling events would be expected on the 
time scale of a typical simulation. 

As can be seen in Figs. II II the asymptotic autocorrela- 
tion time for the stripe order is quite long, approximately 
40 Monte Carlo steps, in the superfluid phase at Kj J = 7 
even for the modest system size L = 16. Critical slowing 
down is expected as the critical superfluid-striped point 
is approached. Although we have not yet attempted to 
extract the corresponding dynamic critical exponent, we 
show, in Fig. results for the Monte Carlo time evolu- 
tion of some relevant critical quantities very close to the 
quantum phase transition. Here the lattice size L = 96 
and the temperature is chosen sufficiently low for ob- 
taining ground state expectation values within statisti- 
cal error. The points shown are averages over "bins" of 
I0 4 Monte Carlo steps, and clearly these bin averages 
are not yet statistically independent; the autocorrelation 
times are several 10 4 Monte Carlo steps. These simu- 
lations did include multi-branch cluster updates. Cur- 
rently, high-precision T — > converged simulations of 
this model close to the superfluid- VBS transition are not 
feasible for L much larger than 100. 
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Note the clear anti-correlations between the stripe or- 
der and superfluid density in Fig.^J| These do not, how- 
ever, give an indication of the order of the phase tran- 
sition between the two phases, as anti-correlations are 
expected at both continuous and first-order transitions. 

IV. DIAGONAL PLAQUETTE INTERACTIONS 

In general, diagonal terms can be added to the 
Hamiltonian Eq. @ without the development of a new 
directed-loop algorithm; only the exit probability tables 
change, due to the modified weights of the relevant di- 
agonal (C) vertices illustrated in Fig. [5] (and symmetry 
related sets). We here consider three different diagonal 
terms: (i) one which enhances or suppresses staggered 
("flippable" by the K-term) plaquettes, (ii) a uniform 
external magnetic (Zeeman) field, and (iii) a staggered 
field. 



A. Flippable-plaquette interaction 

The full spin Hamiltonian including the flippable- 
plaquette interaction is 

(ij) {ijkl) (ijkl) 

where the bond (Bij) and plaquette (Pijkl) operators are 
defined in Eqs. and J5J and Qijki is 1 or for flippable 
and non-flippable plaquettes, respectively, or 

Q m = (1/2 + Sf) (1/2 -S?) (1/2 + (1/2 -Sf) 
+ (1/2 - S?) (1/2 + Sf) (1/2 - SI) (1/2 + Sf) . 

This V interaction is interesting as it produces an exactly 
soluble (Rokhsar-Kivelson^ 3 .) point at J — and —K = 
V, and is similar to the term employed in quantum dimer 
models 3 - (see also Ref.0). Hence, its usefulness making a 
connection between numerical and analytical studies of 
microscopic Hamiltonians is immediately obvious. 

To solve the directed loop equations in the presence of 
a diagonal interaction such as the V term, the general 
procedure is simply to identify, and modify, the relevant 
sets of directed-loop equations which include the vertex 



weighted by V. Here, the diagonal term H\_ a of Eq. i|ll|) 
becomes Hi )0 — Cliju + VQijH- In the SSE formal- 
ism, these diagonal matrix elements are represented as 
C- vertices, each vertex having a weight given by Eq. 115fl . 
The modified vertices of interest here are denoted 

w v = UUT \H a \ 1UT) = (TITI \H a \ TITI) (43) 
= v + c, 

where we have represented the plaquettes of the basis 
state |a) by a list of the spin states in the order ijkl 
corresponding to Fig.^a). The directed-loop equations 
that are relevant to this diagonal term are related to 
those illustrated in Fig. [3J but clearly only the closed 
sets that contain fully-staggered vertices are affected by 
V. The set shown in Fig. [3J does not have any such ver- 
tices and hence the corresponding directed-loop solution 
is the same as with V = 0. In Fig. El we show a closed 
set that is affected by V. Here we have included the 
bounce processes because they can no longer be com- 
pletely excluded. The directed-loop equations for this 
set are written as 

Vn + Vu + V\3 + Wl4 = Wi = C, 

V21 + V 2 2 + V 2 3 + f 24 = W V = V + C, 

V31 + i>32 + v 33 + v 34 = W 3 = J/2, (44) 

U41 + «42 + W43 + «44 = Wi = J/2, 

where Wi , W 3 and W4 are the same as those given before 
in Eq. (|l(j|) . In order to solve these equations, we use the 
same detailed balance requirements, Eq. $T7\. and sym- 
metry arguments as previously to constrain the equations 
and produce a unique solution. Our choice of symmetry 
conditions here correspond to 

Vl2 = «13, 

V22 = v 23l (45) 
W34 = w 44 . 

Two forms of the solutions are needed in order to en- 
sure positive-definite vertex weights for all choices of pa- 
rameters. The first solution is valid for small couplings, 
\V\ < J, and can be formulated without the undesirable 
bounce processes (the right-hand column): 



vii = C - J/2 + V/2, v 12 = J/4 - V/4, vi 3 = J/4 - V/4, v u = 0, 

V21 =C-J/2 + V/2, V22 = J/4 + V/4, v 23 = J/4 + V/4, v 24 = 0, . , 

V31 - J/4 - V/4, W32 - J/4 + V/4, a 33 = 0, v 34 = 0, [W) 

V41 = J/4 - V/4, V42 = J/4 + V/4, a 43 = 0, v 44 = 0. 

where we see that the constant C must be greater than J/2 — V/2 to ensure that Un remains positive. However, 
in order to satisfy the requirement that an is positive for vertex-path sets not affected by the V term, one should 
set C > J/2 for the case V > 0. For V < 0, one must in addition ensure that the weight W v in Eq. l(4*4*|l remains 
positive, requiring a constant C > J/2 + |V|/2. Note that in the limit V — > 0, this equation set is equal to the 
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FIG. 17: A closed set of C and J vertex-paths for which the corresponding directed-loop equations are different for V = and 
V 0. A second set of vertex-paths — the spin-reverse of the one shown here — also gives the same directed-loop solution. 



directed-loop equations previously obtained in Eq. (|18f) . Clearly, when \ V\ > J, a different form of the solutions must 
be constructed in order to ensure that each vertex weight stays positive. In the case V > J > 0, a non-zero bounce 
probability v 2 4 = V — J must be included. A convenient form of the solution in this case is then 

un = C, V12 = 0, V13 = 0, U14 = 0, 

v 2 i = C, v 2 2 = J/2, V23 = J/2, V24, = V - J, , , 

V31 = 0, V32 = J/2, U33 = 0, v u = 0, 

V41 = 0, V42 = J/2, V43 = 0, v M = 0, 

where the bounce process is turned on slowly, i.e. linearly with V — J, ensuring a small bounce probability in the 
algorithm for a moderate range of V larger than the exchange. However, it is clear that the bounce process v 2 4 
becomes negative for negative V (i.e. V < — J < 0) and we hence need a different solution in this case. Again, we 
use a solution that turns the bounce processes on slowly, but in this case two non-zero bounces are required: 



(48) 





= C + J - 


2m 


V12 


-J/2, 


via 


= J/2, 


Vu 


= 2\V\ - 2J, 


''21 


= C + J- 


2|V|, 


V22 


= 0, 


V23 


= 0, 


V2i 


= \v\ - J, 


V31 


= J/2, 




V32 


= 0. 


V33 


= 0, 


V3i 


= 0, 


V41 


= J/2, 




V42 


= 0, 


V43 


= 0, 


V44 


= 0. 



This last equation set imposes the requirement that 
C > 2\V\ - J. In Eqs. gBJ} to gSJl, the actual exit prob- 
abilities for the directed-loop algorithm are obtained in 
the usual way by dividing the matrix elements by the 
vertex weights; P"j — Vij/Wi, where Wi is the relevant 
matrix element. 

Note again that when implementing a diagonal interac- 
tion such as the V term, the only change required in the 
simulation code, relative to the pure J-K model, is the 
probability weights of only the specific relevant vertex- 
paths affected. In the case above for the plaquette V 
term, only the vertex set shown in Fig. 1171 and the re- 
lated set with the other staggered vertex, will use the 
solutions outlined in Eqs. © and ijlB]). All other 

C to J vertex sets which do not contain fully-staggered 
diagonal vertices will use the original solution, Eq. (fT5|> . 



B. Uniform magnetic field 

Perhaps the simplest extension of the J-K Hamilto- 
nian |J3J| is the addition of an external magnetic held. A 



diagonal Zeeman held h coupling to the z-components 
of the S = 1/2 spins is of particular physical impor- 
tance, as it cants the magnetization away from the zero- 
magnetization state, or dopes the system away from half- 
filling in the boson language^ The possibility of decon- 
fined quantum critical points occurring between super- 
fluid and insulating states at commensurate fillings other 
than one- half is also currently a question of interest M> 
The Hamiltonian under study is now 

h = -jJ2 b h - K E p «« - h J2 ( 49 ) 

(ij) (ijkl) 1 

where we restrict h > 0. The diagonal term Hi a of 
Eq. (jl 1|> is modified to include the effects of the field, 

Hi,a = \ {St + S| + S z k +Sf)+ Chjki , (50) 

where this term now produces different matrix elements 
depending on the spins Sf ,Sj,S%, Sf , with an associ- 
ated vertex weight, Eq. JTSJ. We can ensure that each 
weight will remain positive by adjusting C, in particu- 
lar, we write C = h/2 + e, where e > is typically a 



TABLE II: The weight factors for the diagonal vertices in the 
uniform J-K model. 

(Sf S£ SI gf | H a | S'j g| Sj Sf ) weight factor 
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small constant. Representing the relevant plaquette of 
the state |a) by a list of the spin states, we can calcu- 
late the weights of the 16 C- vertices using Eq. (|50|) . The 
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results are summarized in Table ITTI 

In addition, with the Hamiltonian Eq. Ij49(l . the off- 
diagonal plaquette operators H 2 ^ a to Hq, ci in Eq. I|ll|) re- 
main unmodified. In this case, there are now four unique 
sets of directed loop equations for C^C and C^J ver- 
tices that are not related by trivial symmetry operations. 

The first closed set of C and J vertex paths is illus- 
trated in Fig. El with weights ay which now also should 
include the bounce processes an left out in the figure. 
If we recall that the open circle in Fig. [3] denotes a spin 
down, or S z = —1/2, then, the directed-loop equations 
corresponding to this set is now modified from Eq. i|16|) 
to read 

an + a\2 + ai 3 + an = W^ 1 = e, 

021 + a 2 2 + a 23 + «24 = Wi = h/A + e, 

a.3i + a 32 + a 33 + a 34 = W 3 = J/2, (51) 

an + 042 + 043 + «44 = W4 = J/2. 

We use the same detailed balance, Eq. l(TT|l and symmetry 
arguments, Eq. (|45|l . as previously to constrain the equa- 
tion set and produce a unique solution. For fields ft < 4 J, 
we can obtain a solution which contains no bounce pro- 
cesses, 



an = ft/8 - J/2 + e, a 12 = J/4 - ft/16, a 13 

a.21 = h/8 ~ J/2 + e, a 22 = J/4 + ft/16, a 23 

a 3 i = J/4 - ft/16, a 32 = J/4 + ft/16, a 33 

an = J/4 - ft/16, a 42 = J/4 + ft/16, a 43 



J/4 -h/16, ai 4 = 0, 

J/4 + ft/16, a 24 = 0, 

0, a 34 = 0, 

0, 044 = 0, 



(52) 



where, to keep the element an positive-definite, we re- 
quire e > J/2 — ft/8. Again, for ft > 4 J, the form of 
the solutions must change in a non-trivial way in order 
to keep each vertex weight positive. This requirement 
produces a non-zero bounce process, 024, which together 
with the other weights gives the high-field solution 



an 


= e, 


ai2 


-0, 


013 


= 0, 


ai4 


= 0, 


(221 


= e, 


022 


= J/2, 


023 


= J/2, 


024 


= ft/4 


a 3 i 


= 0, 


032 


= J/2, 


«33 


= 0, 


034 


= 0, 


an 


= 0, 


a42 


= J/2, 


043 


= 0, 


044 


= 0. 



(53) 



The second independent closed set of vertex weights 
for the Hamiltonian Eq. (|49|l is obtained by taking the 



spin-reverse of the closed set illustrated in Fig. 03 The 
resulting directed-loop equations then contain the fully 
polarized vertex in Table ^ and are written as 

en + C12 + ci 3 + C14 = Wi = ft + e, 

c 2 i + c 22 + c 23 + c 24 = W^ 1 = 3ft/4 + e, 

C31 + C32 + c 33 + C34 = W 3 = J/2, (54) 

C41 + C42 + C43 + C44 = W 4 = J/2. 
This set is solved in the same way as set a above, employ- 
ing analogous conditions for detailed-balance and vertex 
symmetries. The result is two sets of vertex weights; the 
first, for ft < 4 J is 



en = 7ft/8 - J/2 + e, C12 = J/4 + ft/16, c 13 = J/4 + ft/16, c i4 = 0, 

cai = 7ft/8 - J/2 + e, c 22 = J/4 - ft/16, c 23 = J/4 - ft/16, c 24 = 0, 

C31 = J/4 + ft/16, c 32 = J/4 - ft/16, c 33 - 0, c 34 = 0, 

c 4 i = J/4 + ft/16, c 42 = J/4 - ft/16, c 43 = 0, C44 = 0, 



(55) 
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FIG. 18: A closed set of C and J vertex-paths used in solving the directed-loop equations for the uniform-field J-K-ft model. 



where clearly, to keep Cn > we require e > J/2 — 7h/8. 
We see that this solution is similar in form to the low- 
field solution for vertex set a, Eq. I|52|) . However, the cn 
and C21 terms are modified, and also the other non-zero 
terms have the opposite sign of h (i.e. if a%j = J/4±/i/16, 
then dj = J/4 =p /i/16). The other solution corresponds 
to the large-field case, h > 4 J, and requires the inclusion 
of a non-zero weight for the bounce process c\ 4 , 



cn 


= e + 3/i/4, 


cu 


= J/2, 


Cl3 


= J/2, 


Cl4 


= h/4 


C21 


= e + 3/i/4, 


C22 


= 0, 


C23 


= 0, 


C24 


= 0, 


C31 


= J/2, 


C32 


= 0, 


C33 


= 0, 


C34 


= 0, 


C41 


= J/2, 


C42 


= o, 


C43 


= 0, 


C44 


= 0. 



(56) 

Again, the cn and C21 terms are modified by the different 
vertex sets, and in addition the non-zero bounce process 
has moved to C14 from 024 in Eq. (|53l) . 

The third independent closed set of vertex weights does 
not include the fully spin-up or spin-down matrix ele- 
ments of the previous two, sets a and c. The diagram- 
matic representation is therefore not trivially related to 
these previous cases, and is illustrated in Fig.^J In this 
case, the directed-loop equations are 

dn + du + di3 + d u = W£ = h/4 + e, 

d 2 i + d 22 + d 23 + d 2i = W£ = h/2 + e, 

<*si + d 32 + d 33 + d 3 4 = W 3 = J/2, (57) 

da + d 42 + d 43 + d 44 = W4 = J/2. 

While the detailed balance conditions for this equation 
set is the same as before, Eq. (|17f) . it can be noted that 
the additional symmetry conditions, Eq. (|45l) . do not ap- 
pear in the diagrams here. Although not immediately 
justifiable in terms of symmetry arguments, there is in 
general no reason why the same constrains cannot be 
used to solve equation set d, and therefore we will con- 
tinue to use Eq. (|45|l as it facilitates implementation of 
the algorithm (although there is no guarantee that this 
leads to the most efficient simulation). The low-field so- 
lution (h < 4J) is then given by the same equation set 
as solution a, Eq. (|52|l with all dij — aij except the fol- 
lowing: 

du =d 21 =3/i/8- J/2 + e. (58) 



Here, it is quite obvious that we require e > J/2 — 3h/8 
in order to keep all vertex weights positive. For the high- 
field case (h > 4J), we are forced to have a non-zero 
bounce process, and upon solving we again get an equa- 
tion set similar to solution o, Eq. I|53|l . with the exception 
that 

du = di 2 = e + h/4. (59) 

The final set of vertex-weights used in the uniform- 
field solution is obtained by taking the spin-reverse of the 
closed set illustrated in Fig. 1181 in an analogous manner 
to the way that set c was obtained from set a. The result 
is the directed loop equations given by 

en + era + e 13 + e u = W% = 3h/4 + e, 

e 2 i + e 22 + e 2 3 + e 24 = VF 3 ' 1 = h/2 + e, 

e 3 i + e 32 + e 33 + e 34 = W 3 = J/2, (60) 

e 4 i + e 42 + e 43 + e 44 = W 4 = J/2. 

Again, we employ the detailed balance and symmetry 
conditions discussed above to find a unique low-field so- 
lution for h < 4J, which in this case is the same as the 
solution set c, Eq. I|55|l with the exception that 

en = e 2 i = 5h/8 - J/2 + e, (61) 

where again the bounces have been eliminated, and to get 
en > 0, we need e > J/2-5h/8. The large-field (h > 4J) 
solutions are given by the analogous set, Eq. (|56|l . with 
the exception that 

en = e 2X = e + h/2. (62) 

As before, the above four equations sets, a, c, d and e, 
serve to uniquely define the exit probabilities, given by 
dividing the matrix elements by the vertex weights, e.g. 
Pgj = dij/Wi, where the values of Wi are given either 
by Table |H] above, or by the previously-defined values of 
W 3 — W 4 — J/2. It is important, in the implementation 
of the directed-loop equations, that all diagonal vertices 
are weighted according to their proper equation set. The 
relation of a general C or J vertex to the proper equa- 
tion set in some cases depends on the path that a loop 
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FIG. 19: Sublattice decoration of the two-dimensions square 
lattice, used in constructing the quantum Monte Carlo simu- 
lation of the J-K model with staggered Zeeman field. 



segments takes through the vertex. For example, the 
vertex in 031 is related to ^31, however the path taken 
by the loop in each case results in a C vertex which is 
weighed differently by the field. Also, we note that the 
common element of all of these equation sets is the con- 
stant e. This e must be chosen to keep all of the elements 
an, en, dn, en and their symmetry related weights posi- 
tive definite. The critical condition comes from Eq. (|52[l . 
where an > in all cases for e > J/2 — h/8. It can 
be seen that, if this condition is satisfied, then all of the 
weights ciijdn and en will automatically be positive 
definite, and it is therefore the e that we choose in im- 
plementation of the algorithm. 



C. Staggered magnetic field 

The final set of directed loop solutions that we will 
present in this paper is for the J-K model in a staggered 
Zeeman field. Motivation for this extension of the Hamil- 
tonian comes directly from predictions in the theory of 
deconfined quantum criticalityi^ and its applicability 
to our microscopic model. In short, a staggered Zee- 
man field on our spin model corresponds to a uniform 
Zeeman field that couples to the z component of n in 
the nonlinear-sigma model of relevance. The theory then 
predicts a "split" transition between the VBS and su- 
perfluid phases, with an intermediate phase with neither 
order (but with a "background" field-induced staggered 
magnetization) . 



Hi a of Eq. (|ll(l is modified to include the effects of the 
staggered field, 

H ha = (-1)*^ (S? S; + SI St) + Cl m , (64) 

TABLE III: The weight factors for the diagonal vertices in 
the staggered J-K model. 
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and the other plaquette Hamiltonian terms remain un- 
modified. Keeping C arbitrary for now, we can easily 
calculate the weights for the 16 diagonal (C) vertices. 
The approach we take in constructing the simulation is 
the decorate the lattice with an "A" and "B" sublattice 
in a checkerboard pattern (Fig. ll9|l . The solution to each 
vertex weight in the directed loop equations will have two 
components, one if the vertex happens to fall on an "A" 
plaquette, and another for the same vertex on a "B" pla- 
quette (see Table ITTTfl . 

Turning first to the closed set of C and J diagrams, 
Fig. [S] we construct the directed-loop equations, which 
are now different from the forms Eq. (|16(l and Eq. (|51|l . 



an 
a.21 
a>3i 
an 



ai2 + 013 
CJ22 + 023 

«32 + 033 
042 + (143 



a i4 

«24 
044 



C, 

+h/4 + C, 

J/2, 

J/2. 



(65) 



The modified Hamiltonian is 

H = -jJ^ Bij -K ]T P ijkl -hJ2(-l) Xi+yi S?, (63) 

(<i> (ijkl) i 

where X4 and yi are the Cartesian lattice coordinates of 
the ith spin, and h > 0. The diagonal plaquette term 



Notice that we have suppressed the explicit definition 
used before for the weights, W , in order to simplify no- 
tation. The +~ sign defines the convention that the cor- 
responding term is negative if it falls on an A plaquette, 
and positive if it falls on a B plaquette. We can set the 
bounce processes 014 = a 2 4 = as long as J > h/A, 
giving a solution: 
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FIG. 20: Reference elements of the closed C and J vertex-path diagrams. The element letters (left) and equation set designation 
(right) refer to the relevant vertex weight equation set, Tables llVIIXI 



a n = C- J/2 =f h/8, a 12 = J/4 ± h/16, a 13 = J/4 ± h/16, a u = 0, 

a 21 =C- J/2 T h/8, a 22 = J/4 =F h/16, a 23 = J/4 =F h/16, a 24 = 0, 

a 31 = J/4 ± h/16, a 32 = J/4 =F h/16, a 33 = 0, a 3i = 0, [m) 

a 4i = J/4 ± h/16, a 42 = J/4 =p h/16, a 43 = 0, a 44 = 0. 

Again, the upper symbol of ± or =F refers to the vertex weight on the A sublattice, and the lower symbol corresponds 
to the B sublattice. Note, to keep an > on all plaquettes, we require C > J/2 + h/8. 

Another solution for the A and B sublattices, valid for h > 4J, is obtained by necessarily requiring some non-zero 
bounce probabilities. Consider first the set of equations representing B-plaquettes in Eq. (|65|l . A simple solution can 
be found with a non-zero bounce, a 24 = h/A — J. Note, however, that this sign in the first term in a 24 that makes 
this solution invalid on A-plaquettes (where h —> —h). Thus, for h > 4 J on B plaquettes: 



an 


= c, 


ai2 


= 0, 


ai3 


= 0, 


&14 


= 0, 


a 2 i 


= c, 


a 2 2 


-J/2, 


a 2 3 


= J/2, 


«24 


= h/4 




-o, 


a 32 


-J/2, 




= 0, 


«34 


= 0, 


a 4 i 


= 0, 


042 


= J/2, 


a43 


-o, 


O44 


= 



(67) 

For A-plaquettes, another form of the solution is needed, which is analogous to the solution Eq. (|48|l found for the 
diagonal interaction V < — J < 0. Setting 0,14 = h/2 — 2 J, and a 24 = h/A — J constrains the equations to give the 
solution (h > 4 J on A plaquettes): 



(68) 
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= h/2- 
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a 2 i 


= C + J- 
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= 0. 
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a 24 


= h/4- 


J, 


«31 


= J/2, 




a 32 
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033 


= 0, 


a 34 


= 0, 




0-41 


= J/2, 




a 42 


= 0, 


043 


= 0, 


044 


= 0. 





This imposes the requirement that C > h/2 — J. 

This outlines the basic method of constructing the 
directed-loop probabilities for the staggered magnetic 
field Hamiltonian. The only difficulty in completing the 
procedure is identifying all of the separate sets of closed 
vertex-path diagrams which contribute different vertex 
weights to the directed loop algorithm. Instead of ex- 
plicitly illustrating and solving all of these different sets 
in this case, we simply present the solutions in a more 
concise form. To begin, note that we can abbreviate the 
illustration of the closed sets of C and J vertex-paths if 



we constrain the solutions to obey the same detailed bal- 
ance, Eq. I|17|l . and symmetry arguments, Eq. (|45|) . as 
used throughout this paper. In this case, one only needs 
to know the upper-left (reference) vertex (an in Fig. [SJ) 
in order to uniquely define the entire closed set of C to J 
vertex paths. The rules for constructing the closed set, as 
discussed in section II E, can then be summarized by the 
schematic representation in Fig. [21 which illustrates the 
general relationship between the different vertex weights, 
and their corresponding transformations. 
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TABLE V: Vertex weight equation set II. 



vertex 


A/B (ft < 4J) 


A (ft > 4J) 


B (ft > 4J) 


Sll 


C — J/2 ± ft/8 


C 


C + J — ft/2 


a 12 


J/4 =F ft/ 16 





J/2 


Sl4 








ft/2 - 2J 


S22 


J/4 ± ft/16 


J/2 





S24 





ft/4 - J 


ft/4 ~ J 


S33 











S34 












FIG. 21: Schematic representation of the closed set of C and 
J vertex-paths used in solving the directed-loop equations for 
the J-K model. Blocks with the same shading represent equiv- 
alent vertex weights. 



TABLE IV: Vertex weight equation, set I, for the staggered- 
field J-K model. An example of the starting vertex, an, is 
illustrated in Fig. HJ1 



vertex 


A/B (ft < 4J) 


A (ft > 4J) 


B (ft > 4J) 


an 


C - J/2 =f ft/8 


C+J- ft/2 


C 


ai2 


J/4 ± ft/16 


J/2 





ai4 





ft/2 - 2J 





022 


J/4 =F ft/16 





J/2 


0-24 





ft/4 - J 


ft/4 - J 


033 











034 












We can therefore easily construct an entire closed set of 
C and J vertex paths using Fig. 1211 simply by defining 
the reference vertex. Following this procedure, we see 
that the number of unique vertex probability solution sets 
is narrowed down to six, unrelated by trivial symmetry 
operations. The reference vertices for these six unique 
sets are illustrated in Fig. [2U] The corresponding vertex 
weights are summarized in the equation sets of Tablcs lrvT 
IIXI For example, the first reference vertex of Eq. set I in 
Fig. |2"U1 corresponds to an of Fig. JSJ). The corresponding 
vertex weights appear in Table ITVl and are equivalent to 
the solution sets Eqs. (|(jtifl . (|67l) and |JBSJ. The tables are 
abbreviated to only include unique weights not related 
by the symmetries of Fig. however the full equation 
sets are recovered easily by using this figure or Eqs. i|17|) 
and JJSJl. 



TABLE VI: Vertex weight equation set III. 



vertex A/B (ft < 4 J) A (ft > 4 J) B (ft > 4 J) 

e n C - J/2 ± ft/8 C C + J-h/2 

ei2 J/4 ± ft/16 J/2 

ei4 ft/4 -J ft/4 -J 

e 22 J/4 T ft/16 J/2 

e 24 ft/2-2J 

e 33 

e 34 



TABLE VII: Vertex weight equation set IV. 



vertex A/B (ft < 4J) A (ft > 4J) B (ft > 4J) 

fn C- J/2 =p ft/8 C+J-ft/2 C 

/12 J/4 T- ft/16 J/2 

/14 ft/4 - J ft/4 - J 

/ 22 J/4 ± ft/16 J/2 

/ 24 ft/2-2J 

/ 33 

/ 3 4 



TABLE VIII: Vertex weight equation set V. 



vertex A/B (ft < 4 J) A (ft > 4 J) B (ft > 4 J) 

in C-J/2±3ft/8 C + ft/4 C + J-3ft/4 

ii2 J/4 T ft/1 6 J/2 

ii4 ft/2 - 2J 

i 22 J/4 ± ft/16 J/2 

i 24 ft/4 - J ft/4 - J 

J33 

i 34 



V. DISCUSSION 

In summary, we have developed in this work an ex- 
tensive algorithmic framework for SSE quantum Monte 
Carlo simulations of the S — 1/2 XY model with ring ex- 
change - the J-K model - on a 2D square lattice. In addi- 
tion to outlining the basic representation of the quantum 



TABLE IX: Vertex weight equation set VI. 



vertex A/B (ft < 4 J) A (ft > 4 J) B (ft > 4 J) 

jn C-J/2=F3ft/8 C* + J-3ft/4 C + ft/4 

ji2 J/4 ± ft/16 J/2 

ji4 ft/2-2J 

J22 J/4 T ft/16 J/2 

J24 ft/4 - J ft/4 - J 

J'33 

j 34 00 
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mechanical partition function as a power series expansion 
of plaquette-operators acting on a chosen basis in the 
S z representation, we have developed advanced imple- 
mentations of the directed-loop and multi-branch cluster 
updates, designed to significantly increase algorithm effi- 
ciency in various parameter regimes of the Hamiltonian. 
We have studied the performance of the various updating 
procedures using autocorrelation functions. We have also 
outlined modifications of the directed-loop equations to 
account for extensions of the J-K Hamiltonian to include 
diagonal (potential energy) operators. Although several 
specific Hamiltonian terms are discussed, the procedure 
developed is sufficiently general to allow for easy exten- 
sions to other diagonal interactions. 

TABLE X: Comparison of ground-state energy (per spin) of 
exact diagonalization and SSE quantum Monte Carlo results 
for various parameter values of the J-K Hamiltonian, on a 4x4 
square lattice. The exchange was set to J = 1/2, and simula- 
tions were performed with 50 million Monte Carlo production 
steps each. The staggered field strength is represented by h 3 . 



K/J 


V/J 


hi J 


h s /J 


E eX act 
















-0.5624863 


-0.56249(1) 


1 











-0.6803518 


-0.68034(1) 
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-1.1530991 


-1.15311(2) 
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-0.6239222 


-0.62392(1) 
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-1.2864452 


-1.28643(2) 
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-3 








-0.3983951 


-0.39838(1) 


1 





2 





-0.7434355 


-0.74343(1) 


5 





6 





-1.5000000 


-1.50001(1) 


4 








2 


-1.2547499 


-1.25476(2) 


1 








5 


-1.3859171 


-1.38590(1) 



The last step needed to provide confidence in the rather 
complex implementation of our directed-loop algorithms 
discussed here is to carry out rigorous testing. We do this 
by comparing SSE data with results obtained by exact 
diagonalization of the Hamiltonian. Table Q compares 
ground-state energies obtained in various algorithmic so- 



lution regimes of the quantum Monte Carlo schemes dis- 
cussed here with exact diagonalization results for 4 x 4 
lattices. In the simulations, the temperature Tq was cho- 
sen sufficiently low for there to be no differences, within 
statistical errors, between simulations carried out at To 
and 2Tq. The absence of any detectable differences be- 
tween the exact and SSE results to a relative statistical 
accuracy of w 10~ 5 illustrates the unbiased nature of 
these calculations. 

The SSE algorithm developed here can be extended 
straightforwardly to J-K models with four-spin exchange 
terms on other lattices. For example, implementation 
of the Hamiltonian on the triangular— and kagome^ lat- 
tices is possible for some parameter regimes without be- 
ing hampered by the sign problem. In particular, quan- 
tum Monte Carlo simulations at or near the RK-point 33 
(J = 0, —K = V) are anticipated to make explicit con- 
nection with predictions from analytical theories. Studies 
on these models are underway, and are expected to reveal 
a rich variety of ground state phenomenon. 

In principle, the scheme can be extended also to multi- 
spin interactions on rings with more than four spins, for 
example the XY model with six-spin exchange on the py- 
rochlore lattice— However, it is clear that the directed- 
loop scheme will then become quite complex, and explicit 
solutions of the type we have presented here may not be 
practical. 
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